#MP Moralization and Foreign Policy Attitudes Replication File: Appendix Analyses for Study 1
####Last updated: 2 June 2025

###Analyses carried out using R version 4.4.3 in RStudio version Version XXX on Lenovo ThinkPad X1 Carbon Gen 12 Intel Core Ultra 7 155U, 1700 Mhz, 12 running Windows 11 Enterprize

#### This file replicates material reported in Appendix B. Run "MP_study1_0.R" and "MP_study1_1.R" first for variable coding, functions, packages, and models!


# Appendix Section B.1: Sample Characteristics --------------------------------------------------

round(table(data$man)/sum(table(data$man)), digits=3)

round(table(data$gensilent)/sum(table(data$gensilent)), digits=3)[2]
round(table(data$genboom)/sum(table(data$genboom)), digits=3)[2]
round(table(data$genx)/sum(table(data$genx)), digits=3)[2]
round(table(data$genmill)/sum(table(data$genmill)), digits=3)[2]
round(table(data$genz)/sum(table(data$genz)), digits=3)[2]

round(table(data$university)/sum(table(data$university)), digits=3)[2]

round(table(data$nhwhite)/sum(table(data$nhwhite)), digits=3)[2]

round(table(data$pidreplean)/sum(table(data$pidreplean)), digits=3)[2]
round(table(data$piddemlean)/sum(table(data$piddemlean)), digits=3)[2]
round(table(data$pidtrueind)/sum(table(data$pidtrueind)), digits=3)[2]


round(table(data$regionf, useNA="ifany")/sum(table(data$regionf)), digits=3)[1] #midwest
round(table(data$regionf, useNA="ifany")/sum(table(data$regionf)), digits=3)[2] #northeast
round(table(data$regionf, useNA="ifany")/sum(table(data$regionf)), digits=3)[3] #south
round(table(data$regionf, useNA="ifany")/sum(table(data$regionf)), digits=3)[4] #west

sample.col <- cbind(round(c(table(data$man)/sum(table(data$man)), (table(data$age1)/sum(table(data$age1)))[2], 
                            (table(data$age2)/sum(table(data$age2)))[2], (table(data$age3)/sum(table(data$age3)))[2], 
                            (table(data$age4)/sum(table(data$age4)))[2], (table(data$age5)/sum(table(data$age5)))[2], 
                            (table(data$age6)/sum(table(data$age6)))[2], (table(data$white)/sum(table(data$white)))[2], 
                            (table(data$black)/sum(table(data$black)))[2], (table(data$racenwb)/sum(table(data$racenwb)))[2],
                            (table(data$latinor)/sum(table(data$latinor))), (table(data$regionf)/sum(table(data$regionf)))),digits=3))
target.col <- c(1-0.492, 0.492, 0.119, 0.179, 0.164, 0.160, 0.166, 0.212, 0.754, 0.131, 0.116, 0.837, 0.164, 0.208, 0.171, 0.383, 0.239) #obtained from US Census data

sample.mat <- cbind(sample.col, target.col)
colnames(sample.mat) <- c("Sample", "Population")
rownames(sample.mat) <- c("Not Man", "Man", "18-24", "25-34", "35-44", "45-54", "55-64", "65+", "White", "Black or African-American", "Other race",
                          "Not Hispanic or Latino/a", "Hispanic or Latino/a", "Midwest", "Northeast", "South", "West")

####Table B.1 ###
stargazer(sample.mat)


###References in text:
gentab <- table(data$gender, useNA = "ifany")
sum(head(gentab[-c(1,2)])) # 10 respondents selected either "other" (3) or "prefer not to say" (4)


prog <- table(data$Progress)
sum(head(prog, -1)) # 142 partial responses



# Appendix Section B2: Moralization by Issue Area -------------------------

stargazer(obmod.pooled.1, title="Issue Area Moralization",
          omit.stat=c("LL", "ser", "f","adj.rsq"), style="apsr", digits=2, label="tab:issuemods", se=list(obmod.pooled.1r[,2]),
          omit = c("inc1", "inc2", "inc3", "regionfnortheast", "regionfsouth", "regionfwest", "rel_attendr", "rel_impr"),
          covariate.labels = c("Aid", "Climate", "Cooperation",  "Human Rights", "Security", "Partner",
                               "Adversary", "Man", "Non-Hispanic White", "Gen: Silent", "Gen: Boomer", "Gen: X", "Gen: Millennial",
                               "Republican", "Democrat", "University", "Pol. Interest", "MI", "CI", "Isolationism"), 
          star.cutoffs = c(0.1, 0.05, 0.01), star.char=c("\\dagger", "*", "**"), notes="$^{\\dagger}$p $<$ .1; $^{*}$p $<$ .05; $^{**}$p $<$ .01")




# #### Appendix Section B.3: Issue Salience -------------------------------

###NOTE: Figure B1 (Issue Salience in White House Speeches) generated in Stata; data and .do file available in replication materials ###

##Issue salience ratings

#Calculate means and sds for each policy
hr1salience.mean <- round(mean(data$salhr1, na.rm=TRUE), digits=3)
hr2salience.mean <- round(mean(data$salhr2, na.rm=TRUE), digits=3)
hr1salience.sd <- round(sd(data$salhr1, na.rm=TRUE), digits=3)
hr2salience.sd <- round(sd(data$salhr2, na.rm=TRUE), digits=3)

climate1salience.mean <- round(mean(data$salclimate1, na.rm=TRUE), digits=3)
climate2salience.mean <- round(mean(data$salclimate2, na.rm=TRUE), digits=3)
climate1salience.sd <- round(sd(data$salclimate1, na.rm=TRUE), digits=3)
climate2salience.sd <- round(sd(data$salclimate2, na.rm=TRUE), digits=3)

aid1salience.mean <- round(mean(data$salaid1, na.rm=TRUE), digits=3)
aid2salience.mean <- round(mean(data$salaid2, na.rm=TRUE), digits=3)
aid1salience.sd <- round(sd(data$salaid1, na.rm=TRUE), digits=3)
aid2salience.sd <- round(sd(data$salaid2, na.rm=TRUE), digits=3)

security1salience.mean <- round(mean(data$salsecurity1, na.rm=TRUE), digits=3)
security2salience.mean <- round(mean(data$salsecurity2, na.rm=TRUE), digits=3)
security1salience.sd <- round(sd(data$salsecurity1, na.rm=TRUE), digits=3)
security2salience.sd <- round(sd(data$salsecurity2, na.rm=TRUE), digits=3)

cooperation1salience.mean <- round(mean(data$salcooperation1, na.rm=TRUE), digits=3)
cooperation2salience.mean <- round(mean(data$salcooperation2, na.rm=TRUE), digits=3)
cooperation1salience.sd <- round(sd(data$salcooperation1, na.rm=TRUE), digits=3)
cooperation2salience.sd <- round(sd(data$salcooperation2, na.rm=TRUE), digits=3)

economic1salience.mean <- round(mean(data$saleconomic1, na.rm=TRUE), digits=3)
economic2salience.mean <- round(mean(data$saleconomic2, na.rm=TRUE), digits=3)
economic1salience.sd <- round(sd(data$saleconomic1, na.rm=TRUE), digits=3)
economic2salience.sd <- round(sd(data$saleconomic2, na.rm=TRUE), digits=3)

#Store the information in a matrix
sal.means <- rbind(hr1salience.mean, hr2salience.mean, climate1salience.mean, climate2salience.mean, aid1salience.mean,
                   aid2salience.mean, security1salience.mean,security2salience.mean, cooperation1salience.mean, cooperation2salience.mean,
                   economic1salience.mean, economic2salience.mean)
sal.sds <- rbind(hr1salience.sd, hr2salience.sd, climate1salience.sd, climate2salience.sd, aid1salience.sd,
                 aid2salience.sd, security1salience.sd,security2salience.sd, cooperation1salience.sd, cooperation2salience.sd,
                 economic1salience.sd, economic2salience.sd)
sal.vars <- rbind("Defend Human Rights", "Decent Work", "Fossil Fuels", "Biodiversity", "Digital Inclusion", "Vaccine Programs", 
                  "Counter-terrorism", "Nuclear Sabotage", "United Nations", "International Court of Justice", "Tariffs", "Foreign Direct Investment")

sal.sds <- as.numeric(as.character(sal.sds))
sal.means <- as.numeric(as.character(sal.means))
sal.descript <- cbind(sal.vars, sal.means, sal.sds)
colnames(sal.descript) <- c("Issue", "Mean", "Standard \n Deviation" )

####Table B3: Issue Salience Ratings
xtable(sal.descript, label="tab:issuesal", caption="Issue Salience Ratings") 

#Comparison of mean tests for salience
t.sal.hr <-t.test(data$salhr1, data$salhr2)
t.sal.climate <-t.test(data$salclimate1, data$salclimate2)
t.sal.aid <-t.test(data$salaid1, data$salaid2)
t.sal.coop <-t.test(data$salcooperation1, data$salcooperation2)
t.sal.sec <-t.test(data$salsecurity1, data$salsecurity2)
t.sal.econ <-t.test(data$saleconomic1, data$saleconomic2)

ts.salience <- round(cbind(t.sal.hr$statistic, t.sal.climate$statistic, t.sal.aid$statistic, t.sal.coop$statistic, t.sal.sec$statistic, t.sal.econ$statistic), digits = 3)
ps.salience <- round(cbind(t.sal.hr$p.value, t.sal.climate$p.value, t.sal.aid$p.value, t.sal.coop$p.value, t.sal.sec$p.value, t.sal.econ$p.value), digits = 3)

rbind(ts.salience, ps.salience) #All differences p<0.001

# #### Appendix Section B.4: Subgroup Differences -------------------------


#Men vs women
t.hrgen <- t.test(dat.hrm$dv_moral, dat.hrf$dv_moral) #n.s.
mean.hrm <- mean(dat.hrm$dv_moral, na.rm=TRUE)
mean.hrf <- mean(dat.hrf$dv_moral, na.rm=TRUE)
diff.hrgen <- round(mean.hrm - mean.hrf, digits=3)

t.aidgen <-t.test(dat.aidm$dv_moral, dat.aidf$dv_moral) #p<0.01
mean.aidm <- mean(dat.aidm$dv_moral, na.rm=TRUE)
mean.aidf <- mean(dat.aidf$dv_moral, na.rm=TRUE)
diff.aidgen <- round(mean.aidm - mean.aidf, digits=3)

t.climategen <-t.test(dat.climatem$dv_moral, dat.climatef$dv_moral) #p<0.01
mean.climatem <- mean(dat.climatem$dv_moral, na.rm=TRUE)
mean.climatef <- mean(dat.climatef$dv_moral, na.rm=TRUE)
diff.climategen <- round(mean.climatem - mean.climatef, digits=3)

t.secgen <-t.test(dat.secm$dv_moral, dat.secf$dv_moral) #p<0.01
mean.secm <- mean(dat.secm$dv_moral, na.rm=TRUE)
mean.secf <- mean(dat.secf$dv_moral, na.rm=TRUE)
diff.secgen <- round(mean.secm - mean.secf, digits=3)

t.coopgen <-t.test(dat.coopm$dv_moral, dat.coopf$dv_moral) #p<0.01
mean.coopm <- mean(dat.coopm$dv_moral, na.rm=TRUE)
mean.coopf <- mean(dat.coopf$dv_moral, na.rm=TRUE)
diff.coopgen <- round(mean.coopm - mean.coopf, digits=3)

t.econgen <- t.test(dat.econm$dv_moral, dat.econf$dv_moral) #p<0.01
mean.econm <- mean(dat.econm$dv_moral, na.rm=TRUE)
mean.econf <- mean(dat.econf$dv_moral, na.rm=TRUE)
diff.econgen <- round(mean.econm - mean.econf, digits=3)

#Differences
diffs.gen <- cbind(diff.hrgen, diff.aidgen, diff.climategen, diff.secgen, diff.coopgen, diff.econgen)
ts.gen <- round(cbind(t.hrgen$statistic, t.aidgen$statistic, t.climategen$statistic, t.secgen$statistic, t.coopgen$statistic, t.econgen$statistic), digits = 3)


###Older vs. Younger
t.hrage <- t.test(dat.hrold$dv_moral, dat.hryoung$dv_moral) #n.s.
mean.hrold <- mean(dat.hrold$dv_moral, na.rm=TRUE)
mean.hryoung <- mean(dat.hryoung$dv_moral, na.rm=TRUE)
diff.hrage <- round(mean.hrold - mean.hryoung, digits=3)

t.aidage <-t.test(dat.aidold$dv_moral, dat.aidyoung$dv_moral) #n.s.
mean.aidold <- mean(dat.aidold$dv_moral, na.rm=TRUE)
mean.aidyoung <- mean(dat.aidyoung$dv_moral, na.rm=TRUE)
diff.aidage <- round(mean.aidold - mean.aidyoung, digits=3) 

t.climateage <-t.test(dat.climateold$dv_moral, dat.climateyoung$dv_moral) #n.s.
mean.climateold <- mean(dat.climateold$dv_moral, na.rm=TRUE)
mean.climateyoung <- mean(dat.climateyoung$dv_moral, na.rm=TRUE)
diff.climateage <- round(mean.climateold - mean.climateyoung, digits=3) #-0.004

t.secage <- t.test(dat.secold$dv_moral, dat.secyoung$dv_moral) #p<0.01
mean.secold <- mean(dat.secold$dv_moral, na.rm=TRUE)
mean.secyoung <- mean(dat.secyoung$dv_moral, na.rm=TRUE)
diff.secage <- round(mean.secold - mean.secyoung, digits=3) 

t.coopage <-t.test(dat.coopold$dv_moral, dat.coopyoung$dv_moral) #p<0.01
mean.coopold <- mean(dat.coopold$dv_moral, na.rm=TRUE)
mean.coopyoung <- mean(dat.coopyoung$dv_moral, na.rm=TRUE)
diff.coopage <- round(mean.coopold - mean.coopyoung, digits=3) 

t.econage <- t.test(dat.econold$dv_moral, dat.econyoung$dv_moral) #p<0.01
mean.econold <- mean(dat.econold$dv_moral, na.rm=TRUE)
mean.econyoung <- mean(dat.econyoung$dv_moral, na.rm=TRUE)
diff.econage <- round(mean.econold - mean.econyoung, digits=3) 

#Differences
diffs.age <- cbind(diff.hrage, diff.aidage, diff.climateage, diff.secage, diff.coopage, diff.econage)
ts.age <- round(cbind(t.hrage$statistic, t.aidage$statistic, t.climateage$statistic, t.secage$statistic, t.coopage$statistic, t.econage$statistic), digits = 3)


###Republican vs. Democrat
t.hrparty <- t.test(dat.hrrep$dv_moral, dat.hrdem$dv_moral) #p<0.01
mean.hrrep <- mean(dat.hrrep$dv_moral, na.rm=TRUE)
mean.hrdem <- mean(dat.hrdem$dv_moral, na.rm=TRUE)
diff.hrparty <- round(mean.hrrep - mean.hrdem, digits=3) 

t.aidparty <-t.test(dat.aidrep$dv_moral, dat.aiddem$dv_moral) #p<0.01
mean.aidrep <- mean(dat.aidrep$dv_moral, na.rm=TRUE)
mean.aiddem <- mean(dat.aiddem$dv_moral, na.rm=TRUE)
diff.aidparty <- round(mean.aidrep - mean.aiddem, digits=3) 

t.climateparty <-t.test(dat.climaterep$dv_moral, dat.climatedem$dv_moral) #p<0.01
mean.climaterep <- mean(dat.climaterep$dv_moral, na.rm=TRUE)
mean.climatedem <- mean(dat.climatedem$dv_moral, na.rm=TRUE)
diff.climateparty <- round(mean.climaterep - mean.climatedem, digits=3) 

t.secparty <-t.test(dat.secrep$dv_moral, dat.secdem$dv_moral) #n.s.
mean.secrep <- mean(dat.secrep$dv_moral, na.rm=TRUE)
mean.secdem <- mean(dat.secdem$dv_moral, na.rm=TRUE)
diff.secparty <- round(mean.secrep - mean.secdem, digits=3) 

t.coopparty <-t.test(dat.cooprep$dv_moral, dat.coopdem$dv_moral) #p<0.01
mean.cooprep <- mean(dat.cooprep$dv_moral, na.rm=TRUE)
mean.coopdem <- mean(dat.coopdem$dv_moral, na.rm=TRUE)
diff.coopparty <- round(mean.cooprep - mean.coopdem, digits=3) 

t.econparty <- t.test(dat.econrep$dv_moral, dat.econdem$dv_moral) #n.s.
mean.econrep <- mean(dat.econrep$dv_moral, na.rm=TRUE)
mean.econdem <- mean(dat.econdem$dv_moral, na.rm=TRUE)
diff.econparty <- round(mean.econrep - mean.econdem, digits=3) 

#Differences
diffs.party <- cbind(diff.hrparty, diff.aidparty, diff.climateparty, diff.secparty, diff.coopparty, diff.econparty)
ts.party <- round(cbind(t.hrparty$statistic, t.aidparty$statistic, t.climateparty$statistic, t.secparty$statistic, t.coopparty$statistic, t.econparty$statistic), digits = 3)

###High vs. low Militant Internationalism
thrhawks <- t.test(dat.hrhawk$dv_moral, dat.hrdove$dv_moral) #p<0.01
mean.hrhawk <- mean(dat.hrhawk$dv_moral, na.rm=TRUE)
mean.hrdove <- mean(dat.hrdove$dv_moral, na.rm=TRUE)
diff.hrhawks <- round(mean.hrhawk - mean.hrdove, digits=3)

taidhawks <-t.test(dat.aidhawk$dv_moral, dat.aiddove$dv_moral) #p<0.01
mean.aidhawk <- mean(dat.aidhawk$dv_moral, na.rm=TRUE)
mean.aiddove <- mean(dat.aiddove$dv_moral, na.rm=TRUE)
diff.aidhawks <- round(mean.aidhawk - mean.aiddove, digits=3)

tclimatehawks <-t.test(dat.climatehawk$dv_moral, dat.climatedove$dv_moral) #p<0.01
mean.climatehawk <- mean(dat.climatehawk$dv_moral, na.rm=TRUE)
mean.climatedove <- mean(dat.climatedove$dv_moral, na.rm=TRUE)
diff.climatehawks <- round(mean.climatehawk - mean.climatedove, digits=3)

tsechawks <- t.test(dat.sechawk$dv_moral, dat.secdove$dv_moral) #p<0.01
mean.sechawk <- mean(dat.sechawk$dv_moral, na.rm=TRUE)
mean.secdove <- mean(dat.secdove$dv_moral, na.rm=TRUE)
diff.sechawks <- round(mean.sechawk - mean.secdove, digits=3)

tcoophawks <-t.test(dat.coophawk$dv_moral, dat.coopdove$dv_moral) #p<0.01
mean.coophawk <- mean(dat.coophawk$dv_moral, na.rm=TRUE)
mean.coopdove <- mean(dat.coopdove$dv_moral, na.rm=TRUE)
diff.coophawks <- round(mean.coophawk - mean.coopdove, digits=3)

teconhawks <-t.test(dat.econhawk$dv_moral, dat.econdove$dv_moral) #p<0.01
mean.econhawk <- mean(dat.econhawk$dv_moral, na.rm=TRUE)
mean.econdove <- mean(dat.econdove$dv_moral, na.rm=TRUE)
diff.econhawks <- round(mean.econhawk - mean.econdove, digits=3)

#Differences
diffs.hawk <- cbind(diff.hrhawks, diff.aidhawks, diff.climatehawks, diff.sechawks, diff.coophawks, diff.econhawks)
ts.hawk <- round(cbind(thrhawks$statistic, taidhawks$statistic, tclimatehawks$statistic, tsechawks$statistic, tcoophawks$statistic, teconhawks$statistic), digits = 3)


###High vs. Low Cooperative Internationalism
thrci <- t.test(dat.hrcih$dv_moral, dat.hrcil$dv_moral) #p<0.01
mean.hrcih <- mean(dat.hrcih$dv_moral, na.rm=TRUE)
mean.hrcil <- mean(dat.hrcil$dv_moral, na.rm=TRUE)
diff.hrci <- round(mean.hrcih - mean.hrcil, digits=3)

taidci <- t.test(dat.aidcih$dv_moral, dat.aidcil$dv_moral) #p<0.01
mean.aidcih <- mean(dat.aidcih$dv_moral, na.rm=TRUE)
mean.aidcil <- mean(dat.aidcil$dv_moral, na.rm=TRUE)
diff.aidci <- round(mean.aidcih - mean.aidcil, digits=3)

tclimateci <- t.test(dat.climatecih$dv_moral, dat.climatecil$dv_moral) #p<0.01
mean.climatecih <- mean(dat.climatecih$dv_moral, na.rm=TRUE)
mean.climatecil <- mean(dat.climatecil$dv_moral, na.rm=TRUE)
diff.climateci <- round(mean.climatecih - mean.climatecil, digits=3)

tsecci <- t.test(dat.seccih$dv_moral, dat.seccil$dv_moral) #p<0.01
mean.seccih <- mean(dat.seccih$dv_moral, na.rm=TRUE)
mean.seccil <- mean(dat.seccil$dv_moral, na.rm=TRUE)
diff.secci <- round(mean.seccih - mean.seccil, digits=3)

tcoopci <-  t.test(dat.coopcih$dv_moral, dat.coopcil$dv_moral) #p<0.01
mean.coopcih <- mean(dat.coopcih$dv_moral, na.rm=TRUE)
mean.coopcil <- mean(dat.coopcil$dv_moral, na.rm=TRUE)
diff.coopci <- round(mean.coopcih - mean.coopcil, digits=3)

teconci <- t.test(dat.econcih$dv_moral, dat.econcil$dv_moral) #p<0.01
mean.econcih <- mean(dat.econcih$dv_moral, na.rm=TRUE)
mean.econcil <- mean(dat.econcil$dv_moral, na.rm=TRUE)
diff.econci <- round(mean.econcih - mean.econcil, digits=3)

#Differences
diffs.ci <- cbind(diff.hrci, diff.aidci, diff.climateci, diff.secci, diff.coopci, diff.econci)
ts.ci <- round(cbind(thrci$statistic, taidci$statistic, tclimateci$statistic, tsecci$statistic, tcoopci$statistic, teconci$statistic), digits = 3)


###Iso vs. not ISO
thriso <- t.test(dat.hrisoh$dv_moral, dat.hrisol$dv_moral) #p<0.01
mean.hrisoh <- mean(dat.hrisoh$dv_moral, na.rm=TRUE)
mean.hrisol <- mean(dat.hrisol$dv_moral, na.rm=TRUE)
diff.hriso <- round(mean.hrisoh - mean.hrisol, digits=3)

taidiso <- t.test(dat.aidisoh$dv_moral, dat.aidisol$dv_moral) #p<0.01
mean.aidisoh <- mean(dat.aidisoh$dv_moral, na.rm=TRUE)
mean.aidisol <- mean(dat.aidisol$dv_moral, na.rm=TRUE)
diff.aidiso <- round(mean.aidisoh - mean.aidisol, digits=3)

tclimateiso <- t.test(dat.climateisoh$dv_moral, dat.climateisol$dv_moral) #p<0.01
mean.climateisoh <- mean(dat.climateisoh$dv_moral, na.rm=TRUE)
mean.climateisol <- mean(dat.climateisol$dv_moral, na.rm=TRUE)
diff.climateiso <- round(mean.climateisoh - mean.climateisol, digits=3)

tseciso <- t.test(dat.secisoh$dv_moral, dat.secisol$dv_moral) #n.s.
mean.secisoh <- mean(dat.secisoh$dv_moral, na.rm=TRUE)
mean.secisol <- mean(dat.secisol$dv_moral, na.rm=TRUE)
diff.seciso <- round(mean.secisoh - mean.secisol, digits=3)

tcoopiso <- t.test(dat.coopisoh$dv_moral, dat.coopisol$dv_moral) #p<0.01
mean.coopisoh <- mean(dat.coopisoh$dv_moral, na.rm=TRUE)
mean.coopisol <- mean(dat.coopisol$dv_moral, na.rm=TRUE)
diff.coopiso <- round(mean.coopisoh - mean.coopisol, digits=3)

teconiso <-t.test(dat.econisoh$dv_moral, dat.econisol$dv_moral) #n.s.
mean.econisoh <- mean(dat.econisoh$dv_moral, na.rm=TRUE)
mean.econisol <- mean(dat.econisol$dv_moral, na.rm=TRUE)
diff.econiso <- round(mean.econisoh - mean.econisol, digits=3)


#Differences
diffs.iso <- cbind(diff.hriso, diff.aidiso, diff.climateiso, diff.seciso, diff.coopiso, diff.econiso)
ts.iso <- round(cbind(thriso$statistic, taidiso$statistic, tclimateiso$statistic, tseciso$statistic, tcoopiso$statistic, teconiso$statistic), digits = 3)


#White Non-hispanic vs. Other race or ethnicity
t.hrrace <- t.test(dat.hrnw$dv_moral, dat.hrw$dv_moral) #p<0.01
mean.hrnw <- mean(dat.hrnw$dv_moral, na.rm=TRUE)
mean.hrw <- mean(dat.hrw$dv_moral, na.rm=TRUE)
diff.hrrace <- round(mean.hrnw - mean.hrw, digits=3)

t.aidrace <-t.test(dat.aidnw$dv_moral, dat.aidw$dv_moral) #p<0.01
mean.aidnw <- mean(dat.aidnw$dv_moral, na.rm=TRUE)
mean.aidw <- mean(dat.aidw$dv_moral, na.rm=TRUE)
diff.aidrace <- round(mean.aidnw - mean.aidw, digits=3)

t.climaterace <-t.test(dat.climatenw$dv_moral, dat.climatew$dv_moral) #p<0.01
mean.climatenw <- mean(dat.climatenw$dv_moral, na.rm=TRUE)
mean.climatew <- mean(dat.climatew$dv_moral, na.rm=TRUE)
diff.climaterace <- round(mean.climatenw - mean.climatew, digits=3)

t.secrace <-t.test(dat.secnw$dv_moral, dat.secw$dv_moral) #n.s.
mean.secnw <- mean(dat.secnw$dv_moral, na.rm=TRUE)
mean.secw <- mean(dat.secw$dv_moral, na.rm=TRUE)
diff.secrace <- round(mean.secnw - mean.secw, digits=3)

t.cooprace <-t.test(dat.coopnw$dv_moral, dat.coopw$dv_moral) #p<0.05
mean.coopnw <- mean(dat.coopnw$dv_moral, na.rm=TRUE)
mean.coopw <- mean(dat.coopw$dv_moral, na.rm=TRUE)
diff.cooprace <- round(mean.coopnw - mean.coopw, digits=3)

t.econrace <- t.test(dat.econnw$dv_moral, dat.econw$dv_moral) #p<0.01
mean.econnw <- mean(dat.econnw$dv_moral, na.rm=TRUE)
mean.econw <- mean(dat.econw$dv_moral, na.rm=TRUE)
diff.econrace <- round(mean.econnw - mean.econw, digits=3)

#Differences
diffs.race <- cbind(diff.hrrace, diff.aidrace, diff.climaterace, diff.secrace, diff.cooprace, diff.econrace)
ts.race <- round(cbind(t.hrrace$statistic, t.aidrace$statistic, t.climaterace$statistic, t.secrace$statistic, t.cooprace$statistic, t.econrace$statistic), digits = 3)


### Combine into a single data frame:
diffs.all <- rbind(diffs.gen, ts.gen, diffs.age, ts.age, diffs.party, ts.party, diffs.hawk, ts.hawk, diffs.ci, ts.ci, diffs.iso, ts.iso, diffs.race, ts.race)
diffs.all <- data.frame(diffs.all)
colnames(diffs.all) <- c("Human Rights", "Foreign Aid", "Climate", "Security", "Cooperation", "Economics")
rownames(diffs.all) <- c("Men - Women", "t-statg", "Older - Younger", "t-stats", "Rep - Dem", "t-statp", "High MI - Low MI", "t-stath", "High CI - Low CI", "t-statc", "High Iso - Low Iso", "t-stati", "Non-white - White", "t-statr")


###Table B4: Subgroup Differences in Moral Conviction Across Issue Areas
xtable(diffs.all, label="ttests") #asterisks for p-values <0.05 added manually 


#Statistic referenced in text: Statistically significant correlation between non-white and human rights, climate, economics in OLS models:
rbind(cbind("human rights", round(obmod.hr.2r["nhwhite", "Estimate"], digits=2), round(obmod.hr.2r["nhwhite", "Pr(>|t|)"], digits=3)), #b=-0.03, p<0.01
      cbind("aid", round(obmod.aid.2r["nhwhite", "Estimate"], digits=2), round(obmod.aid.2r["nhwhite", "Pr(>|t|)"], digits=3)), #b=-0.01, p=0.226
      cbind("climate", round(obmod.climate.2r["nhwhite", "Estimate"], digits=2), round(obmod.climate.2r["nhwhite", "Pr(>|t|)"], digits=3)), # b=-0.03, p<0.05
      cbind("cooperation", round(obmod.coop.2r["nhwhite", "Estimate"], digits=2), round(obmod.coop.2r["nhwhite", "Pr(>|t|)"], digits=3)), # b=-0.02, p<0.1
      cbind("security", round(obmod.security.2r["nhwhite", "Estimate"], digits=2), round(obmod.security.2r["nhwhite", "Pr(>|t|)"], digits=3)), #b=-0.01, p=0.519
      cbind("economics", round(obmod.econ.2r["nhwhite", "Estimate"], digits=2), round(obmod.econ.2r["nhwhite", "Pr(>|t|)"], digits=3))) # b=-0.04, p<0.01



# #### Appendix Section B.6: Individual Country Treatment Effects ---------


###Table B5: Effect of Country Treatments, Pooled Model
stargazer(mod.pooled.country, title="Effect of Country Treatments, Pooled Model",
          omit.stat=c("LL", "ser", "f", "adj.rsq"), style="apsr", digits=2, label="tab:country_treats", 
          keep=c("countrytreatCanada", "countrytreatChina", "countrytreatGreat Britain", "countrytreatIran", "countrytreatJapan", "countrytreatRussia", "Constant"), #suppress issue area dummies for space
          se=list(mod.pooled.countryr[,2]), dep.var.labels=c("Moral Conviction"),
          covariate.labels = c("Canada", "China", "Great Britain", "Iran", "Japan", "Russia"), 
          star.cutoffs = c(0.1, 0.05, 0.01), star.char=c("\\dagger", "*", "**"), notes="$^{\\dagger}$p $<$ .1; $^{*}$p $<$ .05; $^{**}$p $<$ .01") #note adds the correct legend; delete the line above it


###Generate data frames with treatment effects and standard errors for each policy model (for Figure B2)

#Human Rights
## Create a data frame with coefficients and standard errors
ateshr <- matrix(NA, nrow=12, ncol=3)
getSE <- function(mod, i){
  b <- mod$coefficients[i, 1]
  se <- mod$coefficients[i,2]
  return(c(b,b-1.96*se,b+1.96*se))
}

ateshr[1,] <- getSE(mod.hr1.countries, 2) #defend human rights
ateshr[2,] <- getSE(mod.hr1.countries, 3)
ateshr[3,] <- getSE(mod.hr1.countries, 4)
ateshr[4,] <- getSE(mod.hr1.countries, 5)
ateshr[5,] <- getSE(mod.hr1.countries, 6)
ateshr[6,] <- getSE(mod.hr1.countries, 7)
ateshr[7,] <- getSE(mod.hr2.countries, 2) #decent work
ateshr[8,] <- getSE(mod.hr2.countries, 3)
ateshr[9,] <- getSE(mod.hr2.countries, 4)
ateshr[10,] <- getSE(mod.hr2.countries, 5)
ateshr[11,] <- getSE(mod.hr2.countries, 6)
ateshr[12,] <- getSE(mod.hr2.countries, 7)

#Climate 
## Create a data frame with coefficients and standard errors
atesclimate <- matrix(NA, nrow=12, ncol=3)
getSE <- function(mod, i){
  #i <- which(names(mod)==var_name)
  b <- mod$coefficients[i, 1]
  se <- mod$coefficients[i,2]
  return(c(b,b-1.96*se,b+1.96*se))
}

atesclimate[1,] <- getSE(mod.climate1.countries, 2) #Fossil Fuels
atesclimate[2,] <- getSE(mod.climate1.countries, 3)
atesclimate[3,] <- getSE(mod.climate1.countries, 4)
atesclimate[4,] <- getSE(mod.climate1.countries, 5)
atesclimate[5,] <- getSE(mod.climate1.countries, 6)
atesclimate[6,] <- getSE(mod.climate1.countries, 7)
atesclimate[7,] <- getSE(mod.climate2.countries, 2) #biodiversity
atesclimate[8,] <- getSE(mod.climate2.countries, 3)
atesclimate[9,] <- getSE(mod.climate2.countries, 4)
atesclimate[10,] <- getSE(mod.climate2.countries, 5)
atesclimate[11,] <- getSE(mod.climate2.countries, 6)
atesclimate[12,] <- getSE(mod.climate2.countries, 7)

#Foreign Aid
## Create a data frame with coefficients and standard errors
atesaid <- matrix(NA, nrow=12, ncol=3)
getSE <- function(mod, i){
  #i <- which(names(mod)==var_name)
  b <- mod$coefficients[i, 1]
  se <- mod$coefficients[i,2]
  return(c(b,b-1.96*se,b+1.96*se))
}

atesaid[1,] <- getSE(mod.aid1.countries, 2) #digital inclusion
atesaid[2,] <- getSE(mod.aid1.countries, 3)
atesaid[3,] <- getSE(mod.aid1.countries, 4)
atesaid[4,] <- getSE(mod.aid1.countries, 5)
atesaid[5,] <- getSE(mod.aid1.countries, 6)
atesaid[6,] <- getSE(mod.aid1.countries, 7)
atesaid[7,] <- getSE(mod.aid2.countries, 2) #vaccine readiness
atesaid[8,] <- getSE(mod.aid2.countries, 3)
atesaid[9,] <- getSE(mod.aid2.countries, 4)
atesaid[10,] <- getSE(mod.aid2.countries, 5)
atesaid[11,] <- getSE(mod.aid2.countries, 6)
atesaid[12,] <- getSE(mod.aid2.countries, 7)


#Security
## Create a data frame with coefficients and standard errors
atessecurity <- matrix(NA, nrow=12, ncol=3)
getSE <- function(mod, i){
  #i <- which(names(mod)==var_name)
  b <- mod$coefficients[i, 1]
  se <- mod$coefficients[i,2]
  return(c(b,b-1.96*se,b+1.96*se))
}

atessecurity[1,] <- getSE(mod.security1.countries, 2) #counter-terrorism
atessecurity[2,] <- getSE(mod.security1.countries, 3)
atessecurity[3,] <- getSE(mod.security1.countries, 4)
atessecurity[4,] <- getSE(mod.security1.countries, 5)
atessecurity[5,] <- getSE(mod.security1.countries, 6)
atessecurity[6,] <- getSE(mod.security1.countries, 7)
atessecurity[7,] <- getSE(mod.security2.countries, 2) #nuclear sabotage
atessecurity[8,] <- getSE(mod.security2.countries, 3)
atessecurity[9,] <- getSE(mod.security2.countries, 4)
atessecurity[10,] <- getSE(mod.security2.countries, 5)
atessecurity[11,] <- getSE(mod.security2.countries, 6)
atessecurity[12,] <- getSE(mod.security2.countries, 7)


#Cooperation
## Create a data frame with coefficients and standard errors
atescoop <- matrix(NA, nrow=12, ncol=3)
getSE <- function(mod, i){
  #i <- which(names(mod)==var_name)
  b <- mod$coefficients[i, 1]
  se <- mod$coefficients[i,2]
  return(c(b,b-1.96*se,b+1.96*se))
}

atescoop[1,] <- getSE(mod.coop1.countries, 2) #international institutions/UN
atescoop[2,] <- getSE(mod.coop1.countries, 3)
atescoop[3,] <- getSE(mod.coop1.countries, 4)
atescoop[4,] <- getSE(mod.coop1.countries, 5)
atescoop[5,] <- getSE(mod.coop1.countries, 6)
atescoop[6,] <- getSE(mod.coop1.countries, 7)
atescoop[7,] <- getSE(mod.coop2.countries, 2) #ICJ
atescoop[8,] <- getSE(mod.coop2.countries, 3)
atescoop[9,] <- getSE(mod.coop2.countries, 4)
atescoop[10,] <- getSE(mod.coop2.countries, 5)
atescoop[11,] <- getSE(mod.coop2.countries, 6)
atescoop[12,] <- getSE(mod.coop2.countries, 7)

#Economic
## Create a data frame with coefficients and standard errors
ateseconomic <- matrix(NA, nrow=12, ncol=3)
getSE <- function(mod, i){
  #i <- which(names(mod)==var_name)
  b <- mod$coefficients[i, 1]
  se <- mod$coefficients[i,2]
  return(c(b,b-1.96*se,b+1.96*se))
}

ateseconomic[1,] <- getSE(mod.economic1.countries, 2) #lifting tariffs
ateseconomic[2,] <- getSE(mod.economic1.countries, 3)
ateseconomic[3,] <- getSE(mod.economic1.countries, 4)
ateseconomic[4,] <- getSE(mod.economic1.countries, 5)
ateseconomic[5,] <- getSE(mod.economic1.countries, 6)
ateseconomic[6,] <- getSE(mod.economic1.countries, 7)
ateseconomic[7,] <- getSE(mod.economic2.countries, 2) #FDI Limits
ateseconomic[8,] <- getSE(mod.economic2.countries, 3)
ateseconomic[9,] <- getSE(mod.economic2.countries, 4)
ateseconomic[10,] <- getSE(mod.economic2.countries, 5)
ateseconomic[11,] <- getSE(mod.economic2.countries, 6)
ateseconomic[12,] <- getSE(mod.economic2.countries, 7)


#Plot each set of treatment effects

#Human Rights
#First, organize and label the data frame
names.coefs <- rbind("Canada", "China", "Great Britain", "Iran", "Japan", "Russia", "Canada", "China", "Great Britain", "Iran", "Japan", "Russia")
dvs.ateshr <- rbind("Defend Human Rights", "Defend Human Rights", "Defend Human Rights", "Defend Human Rights", "Defend Human Rights","Defend Human Rights", "Decent Work", "Decent Work", "Decent Work", "Decent Work", "Decent Work", "Decent Work")
df.ateshr <- data.frame(cbind(names.coefs, ateshr, dvs.ateshr))
colnames(df.ateshr) <- c("Treatment", "Effect", "Low", "High", "Issue")
df.ateshr$Effect <- as.numeric(as.character(df.ateshr$Effect))
df.ateshr$Low <- as.numeric(as.character(df.ateshr$Low))
df.ateshr$High <- as.numeric(as.character(df.ateshr$High))
df.ateshr$Treatment <- factor(df.ateshr$Treatment) 
#Coerce an order for the y-axis
df.ateshr$Treatment <- factor(df.ateshr$Treatment , levels = c("Japan", "Canada", "Great Britain", "Iran","China", "Russia"))

#Plot
p.ateshr <- ggplot(df.ateshr, aes(x=Effect, y=Treatment)) +  
  geom_vline(xintercept = 0, colour="gray70", size=1) + 
  geom_pointrangeh(data=df.ateshr, mapping=(aes(x=Effect, y=Treatment, xmin=Low, xmax=High, linetype=factor(Issue), shape=factor(Issue))), 
                   position = position_dodge2v(height = 0.7), size=I(1.05)) + 
  scale_y_discrete(limits = rev(levels(df.ateshr$Treatment))) + xlim(-0.16, 0.12) +
  scale_linetype_discrete(name="") +
  scale_shape_discrete(name="") +
  ylab("") + xlab("Average Treatment Effect (95% CI)") + ggtitle("Human Rights") +
  coord_cartesian(clip="off") +
  theme_bw() + theme(axis.text.y=element_text(color="black", size=16), axis.text.x=element_text(color="black", size=14), 
                     axis.title.y=element_text(color="black"),
                     axis.title.x=element_text(size=12), plot.title=element_text(color="black", size=14, hjust=0.5),
                     legend.text=element_text(size=14), 
                     legend.position="bottom", plot.margin = unit(c(.5,.5,.5,1.5), "cm")) 

#Climate
#First, organize and label the data frame
dvs.atesclimate <- rbind("Fossil Fuels", "Fossil Fuels", "Fossil Fuels", "Fossil Fuels", "Fossil Fuels","Fossil Fuels", "Biodiversity", "Biodiversity", "Biodiversity", "Biodiversity", "Biodiversity", "Biodiversity")
df.atesclimate <- data.frame(cbind(names.coefs, atesclimate, dvs.atesclimate))
colnames(df.atesclimate) <- c("Treatment", "Effect", "Low", "High", "Issue")
df.atesclimate$Effect <- as.numeric(as.character(df.atesclimate$Effect))
df.atesclimate$Low <- as.numeric(as.character(df.atesclimate$Low))
df.atesclimate$High <- as.numeric(as.character(df.atesclimate$High))
df.atesclimate$Treatment <- factor(df.atesclimate$Treatment) 
#Coerce an order for the y-axis
df.atesclimate$Treatment <- factor(df.atesclimate$Treatment , levels = c("Japan", "Canada", "Great Britain", "Iran","China", "Russia"))

#Plot
p.atesclimate <- ggplot(df.atesclimate, aes(x=Effect, y=Treatment)) +  
  geom_vline(xintercept = 0, colour="gray70", size=1) + 
  geom_pointrangeh(data=df.atesclimate, mapping=(aes(x=Effect, y=Treatment, xmin=Low, xmax=High, linetype=factor(Issue), shape=factor(Issue))), 
                   position = position_dodge2v(height = 0.7), size=I(1.05)) + 
  scale_y_discrete(limits = rev(levels(df.atesclimate$Treatment))) + xlim(-0.16, 0.12) +
  scale_linetype_discrete(name="") +
  scale_shape_discrete(name="") +
  ylab("") + xlab("Average Treatment Effect (95% CI)") + ggtitle("Climate") +
  coord_cartesian(clip="off") +
  theme_bw() + theme(axis.text.y=element_text(color="black", size=16), axis.text.x=element_text(color="black", size=14), 
                     axis.title.y=element_text(color="black"),
                     axis.title.x=element_text(size=12), plot.title=element_text(color="black", size=14, hjust=0.5),
                     legend.text=element_text(size=14), 
                     legend.position="bottom", plot.margin = unit(c(.5,.5,.5,1.5), "cm")) 

#Foreign Aid
#First, organize and label the data frame
dvs.atesaid <- rbind("Digital Inclusion", "Digital Inclusion", "Digital Inclusion", "Digital Inclusion", "Digital Inclusion","Digital Inclusion", "Vaccine Readiness", "Vaccine Readiness", "Vaccine Readiness", "Vaccine Readiness", "Vaccine Readiness", "Vaccine Readiness")
df.atesaid <- data.frame(cbind(names.coefs, atesaid, dvs.atesaid))
colnames(df.atesaid) <- c("Treatment", "Effect", "Low", "High", "Issue")
df.atesaid$Effect <- as.numeric(as.character(df.atesaid$Effect))
df.atesaid$Low <- as.numeric(as.character(df.atesaid$Low))
df.atesaid$High <- as.numeric(as.character(df.atesaid$High))
df.atesaid$Treatment <- factor(df.atesaid$Treatment) 
#Coerce an order for the y-axis
df.atesaid$Treatment <- factor(df.atesaid$Treatment , levels = c("Japan", "Canada", "Great Britain", "Iran","China", "Russia"))

#Plot
p.atesaid <- ggplot(df.atesaid, aes(x=Effect, y=Treatment)) +  
  geom_vline(xintercept = 0, colour="gray70", size=1) + 
  geom_pointrangeh(data=df.atesaid, mapping=(aes(x=Effect, y=Treatment, xmin=Low, xmax=High, linetype=factor(Issue), shape=factor(Issue))), 
                   position = position_dodge2v(height = 0.7), size=I(1.05)) + 
  scale_y_discrete(limits = rev(levels(df.atesaid$Treatment))) + xlim(-0.16, 0.12) +
  scale_linetype_discrete(name="") +
  scale_shape_discrete(name="") +
  ylab("") + xlab("Average Treatment Effect (95% CI)") + ggtitle("Foreign Aid") +
  coord_cartesian(clip="off") +
  theme_bw() + theme(axis.text.y=element_text(color="black", size=16), axis.text.x=element_text(color="black", size=14), 
                     axis.title.y=element_text(color="black"),
                     axis.title.x=element_text(size=12), plot.title=element_text(color="black", size=14, hjust=0.5),
                     legend.text=element_text(size=14), 
                     legend.position="bottom", plot.margin = unit(c(.5,.5,.5,1.5), "cm")) 
p.atesaid

#Security
#First, organize and label the data frame
dvs.atessecurity <- rbind("Counter-terrorism", "Counter-terrorism", "Counter-terrorism", "Counter-terrorism", "Counter-terrorism","Counter-terrorism", "Nuclear Sabotage", "Nuclear Sabotage", "Nuclear Sabotage", "Nuclear Sabotage", "Nuclear Sabotage", "Nuclear Sabotage")
df.atessecurity <- data.frame(cbind(names.coefs, atessecurity, dvs.atessecurity))
colnames(df.atessecurity) <- c("Treatment", "Effect", "Low", "High", "Issue")
df.atessecurity$Effect <- as.numeric(as.character(df.atessecurity$Effect))
df.atessecurity$Low <- as.numeric(as.character(df.atessecurity$Low))
df.atessecurity$High <- as.numeric(as.character(df.atessecurity$High))
df.atessecurity$Treatment <- factor(df.atessecurity$Treatment) 
#Coerce an order for the y-axis
df.atessecurity$Treatment <- factor(df.atessecurity$Treatment , levels = c("Japan", "Canada", "Great Britain", "Iran","China", "Russia"))

#Plot
p.atessecurity <- ggplot(df.atessecurity, aes(x=Effect, y=Treatment)) +  
  geom_vline(xintercept = 0, colour="gray70", size=1) + 
  geom_pointrangeh(data=df.atessecurity, mapping=(aes(x=Effect, y=Treatment, xmin=Low, xmax=High, linetype=factor(Issue), shape=factor(Issue))), 
                   position = position_dodge2v(height = 0.7), size=I(1.05)) + 
  scale_y_discrete(limits = rev(levels(df.atessecurity$Treatment))) + xlim(-0.16, 0.12) +
  scale_linetype_discrete(name="") +
  scale_shape_discrete(name="") +
  ylab("") + xlab("Average Treatment Effect (95% CI)") + ggtitle("Security") +
  coord_cartesian(clip="off") +
  theme_bw() + theme(axis.text.y=element_text(color="black", size=16), axis.text.x=element_text(color="black", size=14), 
                     axis.title.y=element_text(color="black"),
                     axis.title.x=element_text(size=12), plot.title=element_text(color="black", size=14, hjust=0.5),
                     legend.text=element_text(size=14), 
                     legend.position="bottom", plot.margin = unit(c(.5,.5,.5,1.5), "cm")) 

#Cooperation
#First, organize and label the data frame
dvs.atescoop <- rbind("United Nations", "United Nations", "United Nations", "United Nations", "United Nations","United Nations", "International Court", "International Court", "International Court", "International Court", "International Court", "International Court")
df.atescoop <- data.frame(cbind(names.coefs, atescoop, dvs.atescoop))
colnames(df.atescoop) <- c("Treatment", "Effect", "Low", "High", "Issue")
df.atescoop$Effect <- as.numeric(as.character(df.atescoop$Effect))
df.atescoop$Low <- as.numeric(as.character(df.atescoop$Low))
df.atescoop$High <- as.numeric(as.character(df.atescoop$High))
df.atescoop$Treatment <- factor(df.atescoop$Treatment) 
#Coerce an order for the y-axis
df.atescoop$Treatment <- factor(df.atescoop$Treatment , levels = c("Japan", "Canada", "Great Britain", "Iran","China", "Russia"))

#Plot
p.atescoop <- ggplot(df.atescoop, aes(x=Effect, y=Treatment)) +  
  geom_vline(xintercept = 0, colour="gray70", size=1) + 
  geom_pointrangeh(data=df.atescoop, mapping=(aes(x=Effect, y=Treatment, xmin=Low, xmax=High, linetype=factor(Issue), shape=factor(Issue))), 
                   position = position_dodge2v(height = 0.7), size=I(1.05)) + 
  scale_y_discrete(limits = rev(levels(df.atescoop$Treatment))) + xlim(-0.16, 0.12) +
  scale_linetype_discrete(name="") +
  scale_shape_discrete(name="") +
  ylab("") + xlab("Average Treatment Effect (95% CI)") + ggtitle("Institutional Cooperation") +
  coord_cartesian(clip="off") +
  theme_bw() + theme(axis.text.y=element_text(color="black", size=16), axis.text.x=element_text(color="black", size=14), 
                     axis.title.y=element_text(color="black"),
                     axis.title.x=element_text(size=12), plot.title=element_text(color="black", size=14, hjust=0.5),
                     legend.text=element_text(size=14), 
                     legend.position="bottom", plot.margin = unit(c(.5,.5,.5,1.5), "cm")) 
p.atescoop


#Economic
##First, organize and label the data frame
dvs.ateseconomic <- rbind("Tariffs", "Tariffs", "Tariffs", "Tariffs", "Tariffs","Tariffs", "Foreign Investment", "Foreign Investment", "Foreign Investment", "Foreign Investment", "Foreign Investment", "Foreign Investment")
df.ateseconomic <- data.frame(cbind(names.coefs, ateseconomic, dvs.ateseconomic))
colnames(df.ateseconomic) <- c("Treatment", "Effect", "Low", "High", "Issue")
df.ateseconomic$Effect <- as.numeric(as.character(df.ateseconomic$Effect))
df.ateseconomic$Low <- as.numeric(as.character(df.ateseconomic$Low))
df.ateseconomic$High <- as.numeric(as.character(df.ateseconomic$High))
df.ateseconomic$Treatment <- factor(df.ateseconomic$Treatment) 
#Coerce an order for the y-axis
df.ateseconomic$Treatment <- factor(df.ateseconomic$Treatment , levels = c("Japan", "Canada", "Great Britain", "Iran","China", "Russia"))

#Plot
p.ateseconomic <- ggplot(df.ateseconomic, aes(x=Effect, y=Treatment)) +  
  geom_vline(xintercept = 0, colour="gray70", size=1) + 
  geom_pointrangeh(data=df.ateseconomic, mapping=(aes(x=Effect, y=Treatment, xmin=Low, xmax=High, linetype=factor(Issue), shape=factor(Issue))), 
                   position = position_dodge2v(height = 0.7), size=I(1.05)) + 
  scale_y_discrete(limits = rev(levels(df.ateseconomic$Treatment))) + xlim(-0.16, 0.12) +
  scale_linetype_discrete(name="") +
  scale_shape_discrete(name="") +
  ylab("") + xlab("Average Treatment Effect (95% CI)") + ggtitle("Economics") +
  coord_cartesian(clip="off") +
  theme_bw() + theme(axis.text.y=element_text(color="black", size=16), axis.text.x=element_text(color="black", size=14), 
                     axis.title.y=element_text(color="black"),
                     axis.title.x=element_text(size=12), plot.title=element_text(color="black", size=14, hjust=0.5),
                     legend.text=element_text(size=14), 
                     legend.position="bottom", plot.margin = unit(c(.5,.5,.5,1.5), "cm")) 


#####Figure B2: Effect of Seperate Country Treatments on Moral Conviction
dev.new(width=12, height=8)
grid.arrange(p.ateshr, p.atesclimate, p.atesaid, p.atescoop, p.atessecurity, p.ateseconomic, nrow=2)






# ###Appendix Section B.7: Issue Support ----------------------------------

#Combine into a list because stargazer doesn't like long model names:
models.supp <- list(mod.hr.2supp, mod.climate.2supp, mod.aid.2supp, mod.security.2supp, mod.coop.2supp, mod.econ.2supp)

###Table B.6: Salience, Treatments, and Issue Support
stargazer(models.supp, title="Moralization by Salience",
          omit.stat=c("LL", "ser", "f", "adj.rsq"), style="apsr", digits=2, label="tab:issuesupp", 
          se=list(mod.hr.2suppr[,2], mod.climate.2suppr[,2], mod.aid.2suppr[,2],mod.security.2suppr[,2],mod.coop.2suppr[,2],mod.econ.2suppr[,2] ),
          keep = c("hisalience", "treatmentfriend", "treatmentadversary", "pidreplean", "piddemlean", "university", "polint1", "mi1", "ci1", "isolat1", "Constant"), #additional controls suppressed for space
          covariate.labels = c("Higher Salience", "Partner", "Adversary", "Republican", "Democrat", "University", "Pol. Interest",
                               "MI", "CI", "Isolationism"), 
          star.cutoffs = c(0.1, 0.05, 0.01), star.char=c("\\dagger", "*", "**"), notes="$^{\\dagger}$p $<$ .1; $^{*}$p $<$ .05; $^{**}$p $<$ .01") #note replaces the legend for statistical signficance; delete the line above it when formatting table


# ###Appendix Section B.8 -------------------------------------------------

#When accounting for covariates, only observe a gender gap in moral conviction about economics:
rbind(cbind("human rights", round(obmod.hr.2r["man", "Estimate"], digits=2), round(obmod.hr.2r["man", "Pr(>|t|)"], digits=3)), #b=-0.02, p=0.153
      cbind("aid", round(obmod.aid.2r["man", "Estimate"], digits=2), round(obmod.aid.2r["man", "Pr(>|t|)"], digits=3)), #b=0.02, p=0.163
      cbind("climate", round(obmod.climate.2r["man", "Estimate"], digits=2), round(obmod.climate.2r["man", "Pr(>|t|)"], digits=3)), # b=0, p=0.745
      cbind("cooperation", round(obmod.coop.2r["man", "Estimate"], digits=2), round(obmod.coop.2r["man", "Pr(>|t|)"], digits=3)), # b=0.01, p=0.552
      cbind("security", round(obmod.security.2r["man", "Estimate"], digits=2), round(obmod.security.2r["man", "Pr(>|t|)"], digits=3)), #b=0.01, p=0.236
      cbind("economic", round(obmod.econ.2r["man", "Estimate"], digits=2), round(obmod.econ.2r["man", "Pr(>|t|)"], digits=3))) # b=0.03, p<0.01

#Generation gaps discussed in text
#younger generation respondents do not moralize climate more than older generation respondents
round(t.climateage$p.value, digits=2) #p=0.73

#older respondents report more moral conviction about security, cooperation, and economics
cbind(diff.secage, round(t.secage$p.value, digits=3)) # difference = 0.051, p<0.01
cbind(diff.coopage, round(t.coopage$p.value, digits=3)) # difference = 0.027, p<0.01
cbind(diff.econage, round(t.econage$p.value, digits=3)) # difference = 0.029, p<0.01



# ###Appendix Section B.10 ------------------------------------------------

##Calculate means and standard deviations for comparison
mean.ssec <- round(mean(data$socialsecmoral, na.rm=TRUE), digits=3) #social security
mean.stem <- round(mean(data$stemcellmoral, na.rm=TRUE), digits=3) # stem cells

ssecmoral.sd <- round(sd(data$socialsecmoral, na.rm=TRUE), digits=3)
stemmoral.sd <- round(sd(data$stemcellmoral, na.rm=TRUE), digits=3)

#combine
means.dom <- rbind(mean.ssec, mean.stem)
sds.dom <- rbind(ssecmoral.sd, stemmoral.sd)
vars.dom <- rbind("Social Security", "Stem Cells")

sds.dom <- as.numeric(sds.dom)
means.dom  <- as.numeric(means.dom)

moral.dom.descript <- cbind(vars.dom, means.dom, sds.dom)

###Table B8: Moral Conviction about Domestic Issues
xtable(moral.dom.descript, label="tab:domesticmoral", caption="Moral Conviction about Domestic Issues")




# ###Supplemental ---------------------------------------------------------


#Pre-registration: "We will also estimate a pooled model that interacts each issue with the treatment 
# indicators to formally assess evidence that the issue moderates the treatment."

####Issue-level heterogeneity
mod.pooled.3 <- lm(dv_moral ~ treatment*issuearea, data=data_long)
mod.pooled.3r <- coeftest(mod.pooled.3, vcov = vcovCL, cluster=~ResponseId)
mod.pooled.3r 

#Add demographics:
mod.pooled.4 <- lm(dv_moral ~ treatment*issuearea + man + nhwhite + gensilent + genboom + genx + genmill +
                     pidreplean + piddemlean + university + inc1 + inc2 + inc3 + regionf, data=data_long)
mod.pooled.4r <- coeftest(mod.pooled.4, vcov = vcovCL, cluster=~ResponseId)
mod.pooled.4r



